
getwd()
setwd("/Users/Desktop")
setwd("/Users/vicissitude_86/Documents//Abortion/Dataset/spss/spss") #### Mac Pro ###
getwd()

setwd("/Users/Aiden/Documents/cambridge_copy/Replicationcourse/Abortion/Dataset/spss/spss") #### IMAC ####
getwd()

setwd("c:/Users/Aiden/Downloads/Abortion/Abortion/Dataset/spss/spss") #### Desktop ###

library(foreign) # in order to read spss files. 
natsal1 <- read.spss("natsal.por", use.value.labels=TRUE, to.data.frame=TRUE) #.por = SPSS portable File 

################################################################## Table 4 #########################################################
############################################################## LOGISTIC REGRESSION ANALYSIS ########################################

natsal1$manyab1 <- natsal1$MANYAB ##new variable will be added to the dataset. # i need to specify that all the values = 1, so the new value will be NA. 
natsal1$manyab1[natsal1$manyab1==-1] <- NA # So [] selects all the collection of elements of the specified variable. 
natsal1$manyab1
natsal1$manyab1[natsal1$manyab1==99] <- NA
natsal1$manyab1
natsal1$manyab1[natsal1$manyab1>2] <-3
natsal1$manyab1
natsal1$manyab2 <- natsal1$manyab1
table(natsal1$manyab2)
natsal1$manyab2[natsal1$manyab2>1] <- 2 
table(natsal1$manyab2)

natsal2$manyab3 <- natsal1$manyab2
natsal2 <- natsal1[natsal1$manyab2!=0,]
View(natsal2)
table(natsal2$manyab2)
#####################################
####################################################################################################################################

#### Age at time of Interview ####

natsal2[[9]]
abortdage <- natsal2[,c("manyab2", "DAGE")]
table(natsal2$manyab2, natsal2$DAGE)
class(natsal2$DAGE)
#### Age at first sexual experience ####
table(natsal2$AFEXP, natsal2$manyab2)
table(natsal2$AFEXP)
class(natsal2$AFEXP)


#### Index of multiple deprivation #### 

table(natsal2$IMDGBQ, natsal2$manyab2) 

######## Categorise into age at last abortion ######################

natsal2[[1068]]
natsal2$AGEABAL2
table(natsal2$manyab2, natsal2$AGEABAL2)

natsal2$agecat <- natsal2$AGEABAL2
natsal2$agecat[natsal2$agecat == -1] <-NA
natsal2$agecat[natsal2$agecat == 99] <-NA

natsal2$agecat[natsal2$agecat<16] <- "<16" ## with quotations it recognises as text.  

natsal2$agecat[natsal2$agecat == "16" | natsal2$agecat=="17"] <- "16-17"
natsal2$agecat[natsal2$agecat== "18" | natsal2$agecat=="19"] <- "18-19"
natsal2$agecat[natsal2$agecat== "20" | natsal2$agecat=="21"| natsal2$agecat== "22" | natsal2$agecat=="23" | natsal2$agecat== "24"] <- "20-24"
natsal2$agecat[natsal2$agecat== "25" | natsal2$agecat=="26"| natsal2$agecat== "27" | natsal2$agecat=="28" | natsal2$agecat== "29"] <- "25-29"
natsal2$agecat[natsal2$agecat== "30" | natsal2$agecat=="31"| natsal2$agecat== "32" | natsal2$agecat=="33" | natsal2$agecat== "34"] <- "30-24"
natsal2$agecat[natsal2$agecat== "35" | natsal2$agecat=="36"| natsal2$agecat== "37" | natsal2$agecat=="38" | natsal2$agecat== "39"] <- "35-39"
natsal2$agecat[natsal2$agecat>40] <- ">40"

table(natsal2$agecat)
class(natsal2$agecat)
natsal2$agecat <- as.factor(natsal2$agecat)
#### Ethnicity ############################

######################################################################### Again (need to combine categorical variables together)                     
names(natsal1)
natsal1[803]
natsal1[["ETHINIC11"]]           
natsal1$ETHNIC11

table(natsal2$manyab2, natsal2$ETHNIC11)                     

table(natsal2$ETHNIC11)
####################################################################
natsal2$ETHNIC11 <-as.character(natsal2$ETHNIC11) ## this transformed to a character, previously categorial
natsal2$ETHNIC11[natsal2$ETHNIC11 == "pakistani" |natsal2$ETHNIC11=="bangladeshi" ] <-  "Pakistani and Bangladeshi"
natsal2$ETHNIC11[natsal2$ETHNIC11 == "chinese" |natsal2$ETHNIC11=="other asian" | natsal2$ETHNIC11 == "other"] <-"chinese and all others"
natsal2$ETHNIC11 <- as.factor(natsal2$ETHNIC11)
table(natsal2$ETHNIC11)
#### so make sure to change from factor var. to character , it allows us to introduce new levels
### recode var (line 416)
### transform back to a factor , you need to transform back to a factor because character will not be allowed in the model
#####################################################################       



#### Number of natural children ####### 
#####Below count from the paper.                      

natsal1[[41]]                     
table(natsal2$manyab2, natsal2$NOCHILD)
class(natsal2$NOCHILD) #numeric value seems fine.
natsal2$NOCHILD <- as.factor(natsal2$NOCHILD)
#### Education ####

natsal2[[982]]
table(natsal2$manyab2, natsal2$EDUC)
class(natsal2$EDUC)                           
##### Housing Tenure ##########################################

natsal1[[391]]
natsal2$TENURE
table(natsal2$manyab2, natsal2$TENURE)

natsal2$TENURE1 <- natsal2$TENURE
table(natsal2$TENURE1)

natsal2$TENURE1 <-as.character(natsal2$TENURE1) ## this transformed to a character, previously categorial
natsal2$TENURE1[natsal2$TENURE1== "rent from private landlord" |natsal2$TENURE1== "rent from council" | natsal2$TENURE1=="rent from housing association" ] <-  "rent"
natsal2$TENURE1 <- as.factor(natsal2$TENURE1)
table(natsal2$TENURE1) 

#### Contraceptive use at first intercourse ####

natsal1[[1064]]
table(natsal2$manyab2, natsal2$FSPRECG2)
natsal2$FSPRECG2 <- as.factor(natsal2$FSPRECG2)


#### Number of sexual partners####


natsal2[666]
table(natsal2$HETTOT) ##this is still a categorial data, so you have to change it to character before recoding.  
class(natsal2$HETTOT)
natsal2$hettot <- as.character(natsal2$HETTOT)
natsal2$hettot[natsal2$hettot == "0" | natsal2$hettot == "1" | natsal2$hettot == "2" | natsal2$hettot == "3-4" ] <- "< 4" # RECODING VARIABLES, character so you need brackets, even if they look like numbers
natsal2$hettot[natsal2$hettot == "5-9" | natsal2$hettot == "10+"] <- ">5"
natsal2$hettot[natsal2$hettot == "not applicable" | natsal2$hettot == "not answered"] <- NA

table(natsal2$hettot)
natsal2$hettot <- as.factor(natsal2$hettot)
#### Birth prior to first abortion #### 

natsal2[[1002]]
natsal2[[1067]]

natsal2$AGE1CHIL[natsal2$AGE1CHIL == -1] <-NA
natsal2$AGE1CHIL[natsal2$AGE1CHIL == 99] <-NA

natsal2$AGEABALL[natsal2$AGEABALL == -1] <-NA
natsal2$AGEABALL[natsal2$AGEABALL == 99] <-NA

natsal2$birfirab <- natsal2$AGE1CHIL < natsal2$AGEABALL  # (binary category)
table(natsal2$birfirab)
class(natsal2$birfirab)
natsal2$birfirab <- as.factor(natsal2$birfirab)
View(natsal2)

#### MULTIVARIATE REGRESSION ####

natsal3 <- natsal2[,c("ETHNIC11", "TENURE", "EDUC","DAGE", "FSPRECG2","birfirab","hettot","IMDGBQ", "NOCHILD","agecat", "AFEXP", "manyab2") ]   #You can't combine data if it is in character, must change to either factor or numeric. Use class(dataset$variable) to check. 
natsal3$manyab2[natsal3$manyab2 == "1"] <- 0
natsal3$manyab2[natsal3$manyab2 == "2"] <- 1
table(natsal3$manyab2)

logit <- glm(manyab2 ~ 
               ETHNIC11 + TENURE+ EDUC + DAGE + FSPRECG2 + birfirab + hettot + NOCHILD + agecat + AFEXP + IMDGBQ, 
             data=natsal3, 
             family = binomial(link = "logit"))
summary(    logit)
exp(coef(   logit))
exp(confint(logit)) # 95% CI for exponentiated coefficients


################################################ Individual Multivariate Regression ##############################################################################################################################################################################################################################################################

#### Contraceptive use at first intercourse####
FSPRECG2


table(natsal2$FSPRECG2, natsal2$manyab2)
natsal2$FSPRECG2 <- as.factor(natsal2$FSPRECG2)
natsal4 <- natsal2[,c("FSPRECG2","manyab2") ]   #You can't combine data if it is in character, must change to either factor or numeric. Use class(dataset$variable) to check. 
natsal4$manyab2[natsal4$manyab2 == "1"] <- 0
natsal4$manyab2[natsal4$manyab2 == "2"] <- 1
table(natsal4$manyab2)
natsal4$FSPRECG2 <- relevel(natsal4$FSPRECG2, "used condom (poss. with other method(s)")

natsal4$FSPRECG2 <- as.factor(natsal4$FSPRECG2)
table(natsal2$FSPRECG2) 

logit <- glm(manyab2 ~
               FSPRECG2, 
             data=natsal4, 
             family = binomial(link = "logit"))
summary(    logit)
exp(coef(   logit))
exp(confint(logit)) # 95% CI for exponentiated coefficients
hetto
#### Sexual Partners####

table(natsal2$hetto, natsal2$manyab2)
natsal2$hetto <- as.factor(natsal2$hetto)
natsal4 <- natsal2[,c("hetto","manyab2") ]   #You can't combine data if it is in character, must change to either factor or numeric. Use class(dataset$variable) to check. 
natsal4$manyab2[natsal4$manyab2 == "1"] <- 0
natsal4$manyab2[natsal4$manyab2 == "2"] <- 1
table(natsal4$manyab2)
natsal4$hetto <- relevel(natsal4$hetto, "< 4")

table(natsal2$hetto) 

logit <- glm(manyab2 ~
               hetto, 
             data=natsal4, 
             family = binomial(link = "logit"))
summary(    logit)
exp(coef(   logit))
exp(confint(logit)) # 95% CI for exponentiated coefficients
#### Education####

table(natsal2$EDUC, natsal2$manyab2)
natsal4 <- natsal2[,c("EDUC","manyab2") ]   #You can't combine data if it is in character, must change to either factor or numeric. Use class(dataset$variable) to check. 
natsal4$manyab2[natsal4$manyab2 == "1"] <- 0
natsal4$manyab2[natsal4$manyab2 == "2"] <- 1
table(natsal4$manyab2)
natsal4$EDUC <- relevel(natsal4$EDUC, "left school aged 17+")

table(natsal2$EDUC) 

logit <- glm(manyab2 ~
               EDUC, 
             data=natsal4, 
             family = binomial(link = "logit"))
summary(    logit)
exp(coef(   logit))
exp(confint(logit)) # 95% CI for exponentiated coefficients
#### Housing Tenure####

table(natsal2$TENURE1, natsal2$manyab2)
natsal4 <- natsal2[,c("TENURE1","manyab2") ]   #You can't combine data if it is in character, must change to either factor or numeric. Use class(dataset$variable) to check. 
natsal4$manyab2[natsal4$manyab2 == "1"] <- 0
natsal4$manyab2[natsal4$manyab2 == "2"] <- 1
table(natsal4$manyab2)
natsal4$TENURE1 <- relevel(natsal4$TENURE1, "own outright/with mortgage or loan")

table(natsal2$TENURE1) 

logit <- glm(manyab2 ~
               TENURE1, 
             data=natsal4, 
             family = binomial(link = "logit"))
summary(    logit)
exp(coef(   logit))
exp(confint(logit)) # 95% CI for exponentiated coefficients
#### Age of multiple deprivation####

table(natsal2$IMDGBQ, natsal2$manyab2)
natsal4 <- natsal2[,c("IMDGBQ","manyab2") ]   #You can't combine data if it is in character, must change to either factor or numeric. Use class(dataset$variable) to check. 
natsal4$manyab2[natsal4$manyab2 == "1"] <- 0
natsal4$manyab2[natsal4$manyab2 == "2"] <- 1
table(natsal4$manyab2)
natsal4$IMDGBQ <- relevel(natsal4$IMDGBQ, "5th quintile (most deprived)")

table(natsal2$IMDGBQ)

logit <- glm(manyab2 ~
               IMDGBQ, 
             data=natsal4, 
             family = binomial(link = "logit"))
summary(    logit)
exp(coef(   logit))
exp(confint(logit)) # 95% CI for exponentiated coefficients
#### Birth prior to first abortion####

table(natsal2$birfirab, natsal2$manyab2)
natsal4 <- natsal2[,c("birfirab","manyab2") ]   #You can't combine data if it is in character, must change to either factor or numeric. Use class(dataset$variable) to check. 
natsal4$manyab2[natsal4$manyab2 == "1"] <- 0
natsal4$manyab2[natsal4$manyab2 == "2"] <- 1
table(natsal4$manyab2)
natsal4$birfirab <- relevel(natsal4$birfirab, "FALSE")

table(natsal2$birfirab)

logit <- glm(manyab2 ~
               birfirab, 
             data=natsal4, 
             family = binomial(link = "logit"))
summary(    logit)
exp(coef(   logit))
exp(confint(logit)) # 95% CI for exponentiated coefficients
#### Age at last abortion####

table(natsal2$agecat, natsal2$manyab2)
natsal4 <- natsal2[,c("agecat","manyab2") ]   #You can't combine data if it is in character, must change to either factor or numeric. Use class(dataset$variable) to check. 
natsal4$manyab2[natsal4$manyab2 == "1"] <- 0
natsal4$manyab2[natsal4$manyab2 == "2"] <- 1
table(natsal4$manyab2)
natsal4$agecat <- relevel(natsal4$agecat, "20-24")



logit <- glm(manyab2 ~
               agecat, 
             data=natsal4, 
             family = binomial(link = "logit"))
summary(    logit)
exp(coef(   logit))
exp(confint(logit)) # 95% CI for exponentiated coefficients
#### Ethnic group####

table(natsal2$ETHNIC11, natsal2$manyab2)
natsal4 <- natsal2[,c("ETHNIC11","manyab2") ]   #You can't combine data if it is in character, must change to either factor or numeric. Use class(dataset$variable) to check. 
natsal4$manyab2[natsal4$manyab2 == "1"] <- 0
natsal4$manyab2[natsal4$manyab2 == "2"] <- 1
table(natsal4$manyab2)
natsal4$ETHNIC11 <- relevel(natsal4$ETHNIC11, "white") ###### to relevel this category, you must use this. For some reason 



logit <- glm(manyab2 ~
               ETHNIC11, 
             data=natsal4, 
             family = binomial(link = "logit"))
summary(    logit)
exp(coef(   logit))
exp(confint(logit)) # 95% CI for exponentiated coefficients

